This vignette shows an example of applying pathway level to infer cell-type-specific co-expression networks and extracting co-expressed gene modules that are enriched for biological functions in cell types.

1. Load packages and data

library(CSCORE)
library(Seurat)
library(ggraph)
library(gridExtra)
source("../parallelCSCORE/functions_2.R")

In this vignette, we use the single cell RNA-sequencing data on Peripheral blood mononuclear cells (PBMC) from COVID patients and healthy controls from Wilk et al., which were also studied in our manuscript. This data set can be downloaded via the following bash script

#wget https://hosted-matrices-prod.s3-us-west-2.amazonaws.com/Single_cell_atlas_of_peripheral_immune_response_to_SARS_CoV_2_infection-25/blish_covid.seu.rds
cat("dowload")
#> dowload

After downloading blish_covid.seu.rds, we load it into the R session

####Load the gene set and convert into list
geneset<-geneIds(getGmt("../parallelCSCORE/h.all.v2023.2.Hs.symbols.gmt"))
pbmc <- readRDS('../blish_covid.seu.rds')
# Use the original UMI counts stored in Assay 'RNA'
DefaultAssay(pbmc) <- 'RNA'

2. Select cell types and gene sets to study

In this example, we focus on B cells and infer the B cell-specific co-expression network.

pbmc_B = pbmc[,pbmc$cell.type.coarse %in% 'B']

Depending on the biological question of interest, one may choose to study the co-expression network for any gene set. Here, we chose to infer the co-expression network for the genes with meaningful expression levels in B cells (top 5000 among 26361 genes). There are several reasons for our choice:

  1. All genes with moderate to high expression levels provides a comprehensive and unbiased set of genes that could have meaningful biological functions in a cell type.

  2. If the genes have much lower expression levels, it would be statistically more challenging and biologically less interesting to infer their co-expressions, as these genes might have almost all UMI counts equal to 0.

In general, it will be up to the users’s choice to select the gene sets to study. We recommend choosing the gene sets that are of interest to your application.

mean_exp = rowMeans(pbmc_B@assays$RNA@counts/pbmc_B$nCount_RNA)
genes_selected = names(sort.int(mean_exp, decreasing = T))[1:5000]

3. Run Pathway Level to infer cell-type-specific co-expression network on the specified gene set

We further subset the B cells to those from healthy control subjects in order to study B-cell specific co-expression network among healthy control B cells.

metaData<-pbmc_B@meta.data
metaData_Healthy<-metaData[metaData$Status=="Healthy",]
metaData_COVID<-metaData[metaData$Status=="COVID",]
set.seed(666);metaData_COVID<-metaData_COVID[sample(1:dim(metaData_COVID)[1],dim(metaData_Healthy)[1]),]

4. Get the Subset of two groups

pbmc_B.sub<-subset(pbmc_B, cells=c(rownames(metaData_Healthy), rownames(metaData_COVID)))
rm(pbmc_B)

5. Calculate the network, or score of networks

cor_data_score<-calculate_seq_correlation(geneset,pbmc_B.sub,genes_selected,"COVID","Healthy", returnType = "score", scoreMethod = "GSCNCA")
#> [1] "IRLS converged after 5 iterations."
#> [1] "25 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.4981% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1726% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1260% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0758% co-expression estimates were smaller than -1 and were set to -1."
#> I am here line 62
save(cor_data_score, file = "COVIDHD_cor_data_score.RData")

Next, based on the thresholded co-expression matrix, we apply WGCNA to extract co-expressed gene modules. In particular, we use CS-CORE estimates to measure co-expressions for single cell RNA-sequencing data, which replace the Pearson correlations used in traditional WGNCA workflow, that suffer from inflated false positives and attenuation bias on single cell data as demonstrated in our manuscript.

6. Calculate the network, or score of networks for permutations

cor_data_perm_list<-runPermute_seq_correlation(geneset,pbmc_B.sub,genes_selected,"COVID","Healthy", ncores = 1, nPermute=20, returnType = "score", scoreMethod = "GSCNCA")
#> ############################################1700173379   1   perm start
#> covid_557.7029 HIP045.1505 covid_561.474 HIP002.2085 covid_555_2.2338HIP023.707 HIP045.862 covid_555_2.1916 HIP043.3964 covid_557.2732[1] "IRLS converged after 4 iterations."
#> [1] "16 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2287% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0799% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "13 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1825% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0686% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700173763   2   perm start
#> HIP045.1746 covid_558.1097 covid_555_1.1080 HIP043.341 covid_558.1332covid_557.4853 HIP043.642 covid_558.3186 covid_558.825 HIP044.1153[1] "IRLS converged after 5 iterations."
#> [1] "15 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2196% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0842% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1614% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0477% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700174146   3   perm start
#> HIP043.339 HIP044.3911 covid_558.2927 covid_557.2319 HIP023.1427HIP015.459 HIP044.1913 HIP043.2136 HIP043.4410 HIP045.1967[1] "IRLS converged after 4 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1838% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0645% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "15 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2064% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0552% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700174497   4   perm start
#> covid_558.251 covid_555_2.1654 HIP043.833 covid_555_2.2709 covid_557.814HIP023.672 HIP043.99 HIP002.1288 covid_555_1.2234 covid_558.100[1] "IRLS converged after 4 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1761% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0708% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "12 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1892% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0547% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700174876   5   perm start
#> covid_558.2853 HIP045.1983 HIP023.3167 HIP015.366 HIP023.3407HIP015.454 HIP043.5229 HIP023.75 HIP043.3507 covid_558.86[1] "IRLS converged after 4 iterations."
#> [1] "15 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1804% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0454% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2067% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0686% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700175267   6   perm start
#> HIP043.859 covid_558.1461 covid_556.765 covid_555_2.494 HIP044.3906covid_556.1259 covid_557.2808 covid_556.539 covid_561.176 HIP015.583[1] "IRLS converged after 5 iterations."
#> [1] "12 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2268% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1072% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "11 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2228% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0757% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700175618   7   perm start
#> HIP044.424 covid_561.717 HIP044.1280 covid_557.6583 HIP044.3116HIP023.3152 HIP015.1888 HIP044.3486 covid_558.2631 covid_557.5542[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1861% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0527% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1962% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0678% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700175966   8   perm start
#> covid_555_2.400 covid_558.1067 HIP043.5069 covid_558.201 covid_555_1.550covid_558.537 HIP045.1505 HIP045.2468 covid_556.812 HIP045.2181[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1799% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0640% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2262% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0725% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700176310   9   perm start
#> HIP015.278 covid_561.2669 covid_555_1.1166 covid_557.6844 covid_555_2.392HIP043.1829 HIP043.4278 covid_556.1175 covid_558.89 covid_557.6295[1] "IRLS converged after 4 iterations."
#> [1] "8 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2076% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0523% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 6 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2362% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0850% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700176654   10  perm start
#> covid_555_1.1578 HIP044.3396 HIP044.579 covid_555_2.720 covid_555_1.1717covid_555_2.1131 HIP023.2215 covid_558.3106 covid_558.411 HIP043.4254[1] "IRLS converged after 4 iterations."
#> [1] "12 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1675% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0551% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "13 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2304% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0955% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700177001   11  perm start
#> covid_558.3585 covid_556.690 covid_557.6543 covid_555_2.4294 covid_557.1505HIP023.2381 HIP002.924 HIP002.503 covid_556.184 covid_558.2896[1] "IRLS converged after 5 iterations."
#> [1] "11 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2104% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0881% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "7 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1904% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0606% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700177345   12  perm start
#> covid_556.458 covid_556.824 covid_558.682 covid_555_1.912 HIP043.1212covid_556.111 HIP043.5051 HIP023.1947 HIP043.460 covid_558.1870[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2019% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0792% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "8 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2223% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0787% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700177694   13  perm start
#> covid_556.1573 HIP045.499 covid_558.2434 covid_557.654 covid_558.204covid_556.106 covid_556.294 HIP015.2028 HIP044.54 HIP044.2338[1] "IRLS converged after 4 iterations."
#> [1] "17 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1736% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0620% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1863% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0390% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700178067   14  perm start
#> covid_556.1753 covid_558.2460 covid_557.867 HIP044.980 HIP043.2117covid_557.5221 HIP023.2286 covid_555_2.354 HIP043.4490 HIP043.5[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2131% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0848% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "10 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2134% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0835% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700178474   15  perm start
#> HIP023.1466 covid_556.693 covid_557.5408 covid_558.34 HIP044.2121covid_556.82 HIP044.1057 HIP043.4874 HIP023.90 covid_558.326[1] "IRLS converged after 4 iterations."
#> [1] "13 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1932% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0616% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "20 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1858% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0621% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700178852   16  perm start
#> HIP045.100 HIP044.3303 covid_557.5011 HIP023.3207 HIP044.242covid_558.94 covid_558.1448 HIP043.5508 covid_558.3669 HIP044.454[1] "IRLS converged after 4 iterations."
#> [1] "8 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1548% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0439% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "19 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2815% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1374% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700179229   17  perm start
#> covid_555_1.2269 HIP002.258 HIP002.175 covid_561.1511 HIP002.1617HIP044.1033 HIP043.3404 covid_556.423 HIP045.137 covid_558.1295[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2361% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1032% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "13 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2352% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1011% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700179609   18  perm start
#> covid_555_2.1324 covid_561.1393 covid_558.2257 covid_558.3035 HIP044.1467covid_558.3023 HIP043.3778 HIP043.443 covid_558.2720 covid_558.2699[1] "IRLS converged after 5 iterations."
#> [1] "10 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.3014% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1653% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "15 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1728% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0530% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700179987   19  perm start
#> HIP044.2154 HIP043.2906 covid_555_2.2365 covid_556.722 HIP044.1099HIP044.3511 HIP043.5282 HIP043.1236 HIP043.3001 HIP023.3128[1] "IRLS converged after 6 iterations."
#> [1] "13 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2279% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0921% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 5 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2149% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0764% co-expression estimates were smaller than -1 and were set to -1."
#> ############################################1700180350   20  perm start
#> covid_557.4552 covid_558.3414 covid_556.1788 covid_557.2194 covid_555_2.4184HIP015.1161 covid_558.2414 covid_558.1522 covid_555_1.2179 covid_556.184[1] "IRLS converged after 4 iterations."
#> [1] "14 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.2076% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0734% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "16 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1734% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0681% co-expression estimates were smaller than -1 and were set to -1."
save(cor_data_perm_list, file = "COVIDHD_cor_data_hallmark_perm_list_20.RData")

7. Calculate p values

pvalues<-sort(getPvaluesTTest(cor_data_score, cor_data_perm_list))
save(pvalues, file = "pvaluesTtest.RData")

8. Generate the figures for top pathways

gglist <-list()
for(ii in 1:1){
  pathway<-names(pvalues)[ii]  
  networkslist<-network_seq_correlation(geneset[[pathway]],pbmc_B.sub,genes_selected,"COVID","Healthy")
  gglist[[ii]]<-network_plot(networkslist,name1="COVID",name2="Healthy", fileName=paste0("rank_",ii, "_", pathway, "_positivePNH_Healthy_Diff"))
}
#> [1] "IRLS converged after 5 iterations."
#> [1] "25 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.4981% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1726% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1260% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0758% co-expression estimates were smaller than -1 and were set to -1."
#> Warning: package 'dplyr' was built under R version 4.1.3
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:GSEABase':
#> 
#>     intersect, setdiff, union
#> The following object is masked from 'package:graph':
#> 
#>     union
#> The following object is masked from 'package:AnnotationDbi':
#> 
#>     select
#> The following objects are masked from 'package:IRanges':
#> 
#>     collapse, desc, intersect, setdiff, slice, union
#> The following objects are masked from 'package:S4Vectors':
#> 
#>     first, intersect, rename, setdiff, setequal, union
#> The following object is masked from 'package:Biobase':
#> 
#>     combine
#> The following objects are masked from 'package:BiocGenerics':
#> 
#>     combine, intersect, setdiff, union
#> The following object is masked from 'package:gridExtra':
#> 
#>     combine
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
#> Warning: package 'igraph' was built under R version 4.1.3
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:dplyr':
#> 
#>     as_data_frame, groups, union
#> The following object is masked from 'package:GSEABase':
#> 
#>     union
#> The following objects are masked from 'package:graph':
#> 
#>     degree, edges, intersection, union
#> The following object is masked from 'package:IRanges':
#> 
#>     union
#> The following object is masked from 'package:S4Vectors':
#> 
#>     union
#> The following objects are masked from 'package:BiocGenerics':
#> 
#>     normalize, path, union
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union
#> Using "stress" as default layout
#> Using "stress" as default layout
#> Using "stress" as default layout
#> Warning: Using the `size` aesthetic in this geom was deprecated in ggplot2 3.4.0.
#> i Please use `linewidth` in the `default_aes` field and elsewhere instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
for(ii in 2:2){
  pathway<-names(pvalues)[ii]  
  networkslist<-network_seq_correlation(geneset[[pathway]],pbmc_B.sub,genes_selected,"COVID","Healthy")
  gglist[[ii]]<-network_plot(networkslist,name1="COVID",name2="Healthy", fileName=paste0("rank_",ii, "_", pathway, "_positivePNH_Healthy_Diff"))
}
#> [1] "IRLS converged after 5 iterations."
#> [1] "25 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.4981% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1726% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1260% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0758% co-expression estimates were smaller than -1 and were set to -1."
#> Using "stress" as default layout
#> Using "stress" as default layout
#> Using "stress" as default layout
for(ii in 3:3){
  pathway<-names(pvalues)[ii]  
  networkslist<-network_seq_correlation(geneset[[pathway]],pbmc_B.sub,genes_selected,"COVID","Healthy")
  gglist[[ii]]<-network_plot(networkslist,name1="COVID",name2="Healthy", fileName=paste0("rank_",ii, "_", pathway, "_positivePNH_Healthy_Diff"))
}
#> [1] "IRLS converged after 5 iterations."
#> [1] "25 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.4981% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1726% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1260% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0758% co-expression estimates were smaller than -1 and were set to -1."
#> Using "stress" as default layout
#> Using "stress" as default layout
#> Using "stress" as default layout
for(ii in 4:4){
  pathway<-names(pvalues)[ii]  
  networkslist<-network_seq_correlation(geneset[[pathway]],pbmc_B.sub,genes_selected,"COVID","Healthy")
  gglist[[ii]]<-network_plot(networkslist,name1="COVID",name2="Healthy", fileName=paste0("rank_",ii, "_", pathway, "_positivePNH_Healthy_Diff"))
}
#> [1] "IRLS converged after 5 iterations."
#> [1] "25 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.4981% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.1726% co-expression estimates were smaller than -1 and were set to -1."
#> [1] "IRLS converged after 4 iterations."
#> [1] "18 among 5000 genes have negative variance estimates. Their co-expressions with other genes were set to 0."
#> [1] "0.1260% co-expression estimates were greater than 1 and were set to 1."
#> [1] "0.0758% co-expression estimates were smaller than -1 and were set to -1."
#> Using "stress" as default layout
#> Using "stress" as default layout
#> Using "stress" as default layout

9. Show plots

Plot the first pathway.

  gridExtra::grid.arrange(gglist[[1]]$gg1, gglist[[1]]$gg2, gglist[[1]]$ggDiff, ncol=3)

10. Show plots

Plot the second pathway.

  library(ggraph)
  gridExtra::grid.arrange(gglist[[2]]$gg1, gglist[[2]]$gg2, gglist[[2]]$ggDiff, ncol=3)

# 11. Show plots

Plot the third pathway.

  library(ggraph)
  gridExtra::grid.arrange(gglist[[3]]$gg1, gglist[[3]]$gg2, gglist[[3]]$ggDiff, ncol=3)

# 12. Show plots

Plot the fourth pathway.

  library(ggraph)
  gridExtra::grid.arrange(gglist[[4]]$gg1, gglist[[4]]$gg2, gglist[[4]]$ggDiff, ncol=3)

The other Analysis is possible